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ABSTRACT 

We performed N-body + hydrodynamical simulations of the formation and evolution 
• of galaxy groups and clusters in a ACDM cosmology. The simulations invoke star 

■^j- | formation, chemical evolution with non-instantaneous recycling, metal dependent ra- 

■^j- , diative cooling, strong starbursts and (optionally) AGN driven galactic super winds, 

effects of a meta-galactic UV field and thermal conduction. The properties of the 
| galaxy populations in two clusters, one Virgo-like (T~ 3 keV) and one (sub) Coma- 

like (T~6 keV), are discussed. The global star formation rates of the cluster galaxies 
are found to decrease very significantly from redshift z~2 to 0, in agreement with ob- 
servations. The total K-band luminosity of the cluster galaxies correlates tightly with 
fS ^ total cluster mass, and for models without additional AGN feedback, the zero point 

of the relation matches the observed one fairly well. Compared to the observed galaxy 
luminosity function, the simulations nicely match the number of intermediate-mass 
galaxies (— 20< Mb ^—17, smaller galaxies being affected by resolution limits) but they 
show a deficiency of bright galaxies in favour of an overgrown central cD. High resolu- 
tion tests indicate that this deficiency is not simply due to numerical "over-merging" . 

The redshift evolution of the luminosity functions from z=l to is mainly driven 
by luminosity evolution, but also by merging of bright galaxies with the cD. 
■ The colour-magnitude relation of the cluster galaxies matches the observed "red 

& ' sequence" , though with a large scatter, and on average galaxy metallicity increases 

with luminosity. As the brighter galaxies are essentially coeval, the colour-magnitude 
relation results from metallicity rather than age effects, as observed. On the whole, 
a top-heavy IMF appears to be preferably required to reproduce also the observed 
colours and metallicities of the stellar populations. 

Key words: cosmology: theory — cosmology: numerical simulations — galaxies: 
clusters — galaxies: formation — galaxies: evolution 
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1 INTRODUCTION 

Clusters of galaxies are of great interest both as cosmological 
probes and as "laboratories" for studying galaxy formation. 
The mass function and number density of galaxy clusters 
as a function of redshift is a powerful diagnostic for the 
determination of cosmological parameters (see Voit 2004 for 
a recent comprehensive review). Besides, clusters represent 
higher than average concentrations of galaxies, with active 
interaction and exchange of material between them and their 
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environment, as testified by the non-primordial composition 
of the hot surrounding gas (e.g. Matteucci & Vettolani 1988; 
Arnaud et al. 1992; Renzini et al. 1993; Finoguenov et al 
2000, 2001; De Grandi et al. 2004; Baumgartner et al. 2005). 

There are observational and theoretical arguments in- 
dicating that clusters are not fair samples of the average 
global properties of the Universe: the morphological mix- 
ture of galaxies in clusters is significantly skewed toward 
earlier types with respect to the field population, implying 
star formation histories peaking at higher redshifts than is 
typical in the field (Dressier 1980; Goto et al. 2003; Kodama 
& Bower 2001; see also Section 3); this is qualitatively in 
line with the expectation that high density regions such as 
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clusters, in a hierarchical bottom-up cosmological scenario, 
evolve at an "accelerated" pace with respect to the rest of 
the Universe (Bower 1991; Diaferio et al. 2001; Benson et al. 
2001). Although clusters are somewhat biased sites of galaxy 
formation, they present the advantage of being bound struc- 
tures with deep potential wells, likely to retain all the mat- 
ter that falls within their gravitational influence; henceforth, 
they represent well-defined, self-contained "pools", where 
one can aim at keeping full track of the process of galaxy 
formation and evolution, and of the global interplay between 
galaxies and their environment. 

The physics of clusters of galaxies has thus received 
increasing attention in the past decade, benefiting from a 
number of X-ray missions measuring the emission of the 
hot intra-cluster gas (e.g., ASCA, ROSAT, XMM, Chandra) 
as well as from extensive optical/NIR surveys probing the 
distribution of galaxies and their properties (e.g. MORPHS, 
SDSS, 2MASS). Understanding cluster physics is also crucial 
to reconstruct, from the observed X-ray luminosity function 
and temperature distribution, the intrinsic mass function of 
clusters as a function of redshift, which is a quantity of pro- 
found cosmological interest (Voit 2004). 

The baryonic mass in clusters is largely in the form of 
a hot intra-cluster medium (ICM), which dominates by a 
factor of 5-10 over the stellar mass (Arnaud et al. 1992; 
Lin, Mohr & Stanford 2003). Consequently, early theoret- 
ical work and numerical simulations concentrated on pure 
gas dynamics when modelling clusters. Recently however, 
attention has focused also on the role of galaxy and star 
formation, and related effects. Star formation locks-up low 
entropy gas, and supplies thermal and kinetic energy to the 
surrounding medium via supernova explosion and shell ex- 
pansion; both processes likely contribute to the observed 
"entropy floor" in low-mass clusters, and the corresponding 
breaking of the scaling relations expected from pure grav- 
itational collapse physics (Voit et al. 2003). Besides, star 
formation is accompanied by the production of new metals 
and the chemical enrichment of the environment; the con- 
siderable (about 1/3 solar) metallicity of the ICM indicates 
that a significant fraction of the metals produced — compa- 
rable or even larger than the fraction remaining within the 
galaxies — is dispersed into the intergalactic medium (Ren- 
zini 2004), affecting the cooling rates of the intra-cluster 
gas. It is thus clear that the hydrodynamical evolution of 
the hot ICM is intimately connected to the formation and 
evolution of cluster galaxies. Only recently, however, due 
to advances in computing capabilities as well as in detailed 
physical modelling, cluster simulations have reached a level 
of sophistication adequate to trace star formation and re- 
lated effects in individual galaxies, and the chemical enrich- 
ment of the ICM by galactic winds, in a reasonably realistic 
way (Valdarnini 2003; Tornatore et al. 2004). 

Indeed so far theoretical predictions of the properties of 
cluster galaxy populations within a fully cosmological con- 
text, have been mainly derived by means of semi-analytical 
models. High resolution, purely N-body cosmological simu- 
lations of the evolution of the collisionless dark matter com- 
ponent, are combined with semi-analytical recipes describ- 
ing galaxy formation and related physics (such as chemical 
enrichment, stellar feedback and exchange of gas and metals 
between galaxies and their environment; see e.g. De Lucia, 
Kauffmann & White 2004 for a recent reference); with such 



schemes, the evolution of galaxies is "painted" on top of 
that of the simulated dark matter haloes and sub-haloes. 
The advantage of this technique is that very high resolution 
can be attained, since pure N-body simulations can han- 
dle larger numbers of particles than hydrodynamical simu- 
lations; besides, a wide range of parameters can be explored 
for baryonic physics (e.g. star formation and feedback effi- 
ciency, Initial Mass Function, etc.). 

In this paper, we present for the first time (to our knowl- 
edge) an analysis of the properties of the galaxy population 
of clusters as predicted directly from cosmological simula- 
tions including detailed baryonic physics, gas dynamics and 
galaxy formation and evolution. The resolution for N-body 
+ hydrodynamical simulations cannot reach the level of the 
purely N-body simulations that constitute the backbone of 
semi-analytical models, so we cannot resolve galaxies at the 
faint end of the luminosity function (Mb>-16). On the other 
hand, our simulations have the advantage of describing the 
actual hydrodynamical response of the ICM to star forma- 
tion, stellar feedback and chemical enrichment. Although 
some uncertain physical processes still necessarily rely on pa- 
rameters (like the star formation efficiency and the feedback 
strength, see Section 2), once these are chosen, the interplay 
between cluster galaxies and their environment follows in a 
realistic fashion, as part of the global cosmological evolution 
of the cluster. 

Moreover, the intimate relation between stellar Initial 
Mass Function (IMF), stellar luminosity, chemical enrich- 
ment, supernova energy input, returned gas fraction and gas 
flows out of/into the galaxies is included self-consistently in 
the simulations (while sometimes these are treated as ad- 
justable, independent parameters in semi-analytical mod- 
els). The properties we obtain for the galaxy populations 
(notably, global star formation rates, luminosity functions 
and colour-magnitude relations) are then the end result of 
ab-initio simulations, with a far minor degree of parameter 
calibration than in semi-analytical schemes. 

In a standard ACDM cosmology, we have performed 
N-body + hydrodynamical (SPH) simulations of the forma- 
tion and evolution of clusters of different mass, on scales of 
groups to moderately rich clusters (emission-weighted tem- 
perature from 1 to 6 keV). In Paper I of this series (Romeo 
et al. 2004, in preparation) we analyze the properties and 
distribution of the hot ICM in the simulated clusters, and 
discuss the effects that star formation and related baryonic 
physics have, on the predicted X-ray properties of the hot 
gas. Several sets of simulations have been carried out, assum- 
ing different IMFs and feedback prescriptions (see Paper I 
and Section 2). The chemical and X-ray properties of the 
ICM are best reproduced assuming a fairly top-heavy IMF 
and a high, though not extreme, feedback (super-wind) ef- 
ficiency (the simulations marked, hereafter, as AY-SW; see 
Paper I and Table 1). 

In this Paper II we focus instead on the properties of 
the galaxy population in our simulated richer clusters (with 
temperatures between 3 and 6 keV), where the number of 
identified galaxies is statistically significant. We will mainly 
discuss the results from the AY-SW simulations, favoured by 
the resulting properties of the ICM; results from simulations 
with different input physics (IMF, wind efficiency, preheat- 
ing) are also discussed for comparison, where relevant. 

In Section 2 we briefly introduce the code and the sim- 
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ulations (full details are given in Paper I), as well as the 
procedure to identify cluster galaxies in the simulations. In 
Section 3 we discuss the global star formation histories of 
cluster galaxies, and in Section 4 we determine global lumi- 
nosities of the clusters simulated with different prescriptions 
for baryonic physics. In Section 5 and 6 we discuss luminos- 
ity functions and colour-magnitude relations of the galaxy 
population in our clusters, and, finally, in Section 7 we sum- 
marize our results. 



2 THE SIMULATIONS 

The code used for the simulations is a significantly improved 
version of the TreeSPH code we used previously for galaxy 
formation simulations (Sommer-Larsen, Gotz & Portinari 
2003). Full details on the code and the simulations are given 
in Paper I, here we recall the main upgrades over the pre- 
vious version. (1) In lower resolution regions an improve- 
ment in the numerical accuracy of the integration of the 
basic equations is obtained by incorporating the "conserva- 
tive" entropy equation solving scheme of Springel & Hern- 
quist (2002). (2) Cold high-density gas is turned into stars 
in a probabilistic way as described in Sommer-Larsen et 
al. (2003). In a star- formation event a SPH particle is con- 
verted fully into a star particle. Non-instantaneous recycling 
of gas and heavy elements is described through probabilistic 
"decay" of star particles back to SPH particles as discussed 
by Lia et al. (2002a). In a decay event a star particle is 
converted fully into a SPH particle. (3) Non-instantaneous 
chemical evolution tracing 10 elements (H, He, C, N, O, 
Mg, Si, S, Ca and Fe) has been incorporated in the code 
following Lia et al. (2002a, b); the algorithm includes su- 
pernova? of type II and type la, and mass loss from stars 
of all masses. For the simulations presented in this paper, 
a Salpeter (1955) IMF and an Arimoto-Yoshii (1987) IMF 
were adopted, both with with mass limits [0.1-100] M®. (4) 
Atomic radiative cooling is implemented, depending both 
on the metallicity of the gas (Sutherland & Dopita 1993) 
and on the meta-galactic UV field, modelled after Haardt 
& Madau (1996). (5) Star-burst driven, galactic super-winds 
are incorporated in the simulations. This is required to expel 
metals from the galaxies and reproduce the observed levels of 
chemical enrichment of the ICM. A burst of star formation is 
modelled in the same way as the "early bursts" of Sommer- 
Larsen et al. (2003): the energy released by SNII explosions 
goes initially into the interstellar medium as thermal energy, 
and gas cooling is locally halted to reproduce the adiabatic 
super-shell expansion phase; a fraction of the supplied en- 
ergy is subsequently converted (by the hydro code itself) 
into kinetic energy of the resulting expanding super-winds 
and/or shells. The super-shell expansion also drives the dis- 
persion of the metals produced by type II supernovae (while 
metals produced on longer timescales are restituted to the 
gaseous phase by the "decay" of the corresponding star par- 
ticles, see point 2 above). The strength of the super- winds is 
modelled via a free parameter / w i n( j which determines how 
large a fraction of the new-born stars partake in such burst- 
ing, super-wind driving star formation. We find that in order 
to get an iron abundance in the ICM comparable to observa- 
tions, / w ind > 0.5 and at a top-heavy IMF must be used. (6) 
Thermal conduction was implemented in the code following 



Cleary & Monaghan (1999), with the addition that effects of 
saturation (Cowie & McKee 1977) were taken into account. 

The groups and clusters were drawn and re-simulated 
from a dark matter (DM)-only cosmological simulation run 
with the code FLY (Antonuccio et al, 2003), for a stan- 
dard flat A Cold Dark Matter cosmological model (h — 0.7, 
fio = 0.3, erg = 0.9) with 150ft" 1 Mpc box-length. When 
re-simulating with the hydro-code, baryonic particles were 
"added" to the original DM ones, which were split according 
to a chosen baryon fraction fb = 0.12. This results in particle 
masses of m gas =m*= 2.5 ■ 10 8 and mDM=l-8 ■ 10 9 ft _1 MQ. 
Gravitational (spline) softening lengths of e gas =e,=2.8 and 
e DM =5.4 h~ kpc, respectively, were adopted. The gravity 
softening lengths were fixed in physical coordinates from 
z—6 to 2=0 and in comoving coordinates at earlier times. 
Particle numbers are in the range 4-10 5 -10 6 SPH+DM par- 
ticles. To test for numerical resolution effects one simulation 
was run with eight times higher mass and two times higher 
force resolution, with 2.3 ■ 10 6 SPH+DM particles having 
m gas =m*= 3.1 • 10 7 , mDM=2.3 ■ 10 8 /i _1 M Q and softening 
lengths of 1.4, 1.4 and 2.7 /i -1 kpc, respectively fSection l5.1ll . 

We selected from the cosmological simulation some 
fairly relaxed systems, not undergoing any major mergers 
since z < 1: two galaxy groups, two small mass clusters, one 
"Virgo-like" cluster (T ~ 3 keV) and a "Mini-Coma" (T ~ 6 
keV) cluster. TreeSPH re-simulations were run with different 
super-wind (SW) prescriptions, IMFs, either Salpeter (Sal) 
or Arimoto-Yoshii (AY), with or without thermal conduc- 
tion, with or without preheating and, as mentioned above, 
at two different numerical resolutions. 

As our reference runs we shall consider simulations with 
an AY IMF, with 70% of the energy feedback from super- 
novae type II going into driving galactic super-winds, zero 
conductivity, and at normal resolution; these will be denoted 
"AY-SW". Some other simulations were run with galac- 
tic super-wind feedback two times (SWx2) and four times 
(SWx4) as energetic as is available from supernovae; the ad- 
ditional amount of energy is assumed to come from AGN 
activity. Others were of the AY-SW type, either with addi- 
tional preheating at z~3 (preh.), as discussed by Tornatore 
et al. (2003), or with thermal conduction included (COND), 
assuming a conductivity of 1/3 of the Spitzer value (e.g., 
Jubelgas, Springel & Dolag 2004). Finally, one series of sim- 
ulations was run with a Salpeter IMF and only early (z>4), 
strong feedback, as in the galaxy formation simulations of 
Sommer-Larsen et al. (2003); as this results in overall fairly 
weak feedback we denote this as: Sal-WFB. We refer the 
reader to Papers I and III, and also Sommer-Larsen (2004), 
for more details. 

The general features and main results for all the simula- 
tions are presented in the companion Paper I. In the present 
paper we mainly discuss the properties of the galaxy pop- 
ulations in the simulated "Virgo" and "Mini-Coma" clus- 
ters, which are in Paper III referred to as clusters "CI" and 
"C2" , respectively. We focus on these two largest simulated 
objects since they contain a significant number of identified 
individual galaxies. The prescriptions adopted for the re- 
simulations of these two clusters are summarized in Table 1. 
Lower mass objects are only included in the analysis of the 
luminosity scaling in Section 4. 

We consider the AY-SW simulation as our "standard" 
run. These simulations provide a satisfactory overall match 
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to the observed properties of the ICM (see Paper I, and also 
Sommer-Larsen 2004). For "Virgo" we also have simulations 
with higher feedback efficiency, with pre-heating or with the 
Salpeter IMF (weak or normal feedback); the corresponding 
results for the galaxy population are shown for comparison, 
where relevant (e.g., Figs. Uland llO^ . 

Because of its high computational demand, the "Coma" 
system was run only with the standard AY-SW prescription, 
with (COND) or without thermal conduction. 

2.1 The identification of individual galaxies 

The stellar contents of both clusters are characterized by a 
massive, central dominant (cD) elliptical galaxy surrounded 
by other galaxies orbiting in the main cluster potential, and 
embedded in an extended envelope of tidally stripped intra- 
cluster stars, unbound from the galaxies themselves. Obser- 
vationally, the brightest galaxy in clusters is usually located 
at or near the cluster centre, where other (massive) galaxies 
likely tend to sink and merge by dynamical friction, mainly 
against the dark matter. Our simulated cDs appear very 
dominant, containing slightly more than half of all the star 
particles in the cluster, in this being more prominent than 
the cDs observed in real clusters. This is related to a "central 
cooling" problem common in models and hydrodynamical 
simulations, both SPH and Adaptive Mesh Refinement (e.g. 
Ciotti et al. 1991; Fabian 1994; Knight & Ponman 1997; 
Suginohara & Ostriker 1998; Lewis et al. 2000; Tornatore 
et al. 2003; Nagai & Kravtsov 2004): after the cluster has 
recovered from the last major merging events, a steady cool- 
ing flow is established at the center of the cluster, being only 
partially attenuated by the strong feedback. The cooled-out 
gas is turned into stars at the base of the cooling flow (r < 10 
kpc). This fairly young stellar population, accumulating at 
the centre of the cD, contributes very significantly to the 
total luminosity, and increases the total stellar mass of the 
cD by of the order a factor of two. In semi-analytical mod- 
els, one can alleviate this problem by artificially quenching 
radiative gas cooling in galactic haloes more massive than 
350 km/sec (Kauffmann et al. 1999). However, over-cooling 
in cluster cores is not merely a technical problem in numer- 
ical simulations, but a problem with the physics in the cen- 
tral cluster regions. XMM and Chandra observations have 
revealed that cooling occurs only down to temperatures of 
about 1-2 keV, so that the former "cooling flow" regions are 
now known to be instead "cool cores" (e.g. Molendi & Piz- 
zolato 2001; Peterson et al. 2001, 2003; Tamura et al. 2001; 
Matsushita et al. 2002) ; the mechanisms that prevent the gas 
from cooling further and finally form stars at a high rate, 
are presently not clearly identified (central AGN feedback, 
magnetohydrodynamic effects, thermal conduction and/or 
heating through gravitational interactions are some candi- 
dates; e.g. Ciotti & Ostriker 2001; Narayan & Medvedev 
2001; Makishima et al. 2001; Fabian, Voigt & Morris 2002; 
Churazov et al. 2002; Kaiser & Binney 2003; Fujita, Suzuki 
& Wada 2004; Cen 2005). 

The cD and its envelope of intracluster stars formed in 
our simulations are discussed in the companion Paper III 
of this series (Sommer-Larsen, Romeo & Portinari, 2005). 
Here, we study the properties of the rest of the galaxy pop- 
ulation in the simulated clusters. The galaxies are identified 
in the simulations by means of the procedure detailed in 



Paper III. Visual inspection of the z=0 frames shows that 
the stars in all galaxies (except the cD) are typically lo- 
cated within 10-15 kpc from the respective galactic centers. 
In a cubic grid of cube-length 10 kpc, we identify all cells 
containing at least 2 star particles. Each of these is then em- 
bedded in a larger cube of length 30 kpc; if this larger cube 
contains a minimum of N m i n =7 gravitationally bound star 
particles, the system is labelled as a potential galaxy. Finally, 
among the various potential galaxies effectively identifying 
the same system, we classify as the galaxy the one contain- 
ing the largest number of star particles. With the mass res- 
olution of the simulations, N m in corresponds to a stellar 
mass of 2.5 x 10 9 M@ and an absolute B-band magnitude 
of Mb ~ —16, of the order of that of the Large Magellanic 
Cloud. The galaxy identification algorithm is adequately ro- 
bust, as long as A r m i„~7-10 star particles. Note though that 
galaxies will contain gas and dark matter particles as well, 
in the case of small galaxies often considerably more dark 
matter particles, because most of the baryons have been ex- 
pelled by galactic super-winds. 

Using this galaxy identification procedure, we identify 
for the standard runs 42 galaxies in "Virgo" and 212 in 
"Coma" , respectively. 



2.2 The computation of the luminosity 

We assign luminosities and colours to the galaxies identified 
in our simulations, as the sum of the luminosities of the 
relevant star particles in the various passbands. 

Each star particle represents a Single Stellar Popula- 
tion (SSP) of total mass of 3.6 x 10 8 Mq, with individual 
stellar masses distributed according to a particular IMF (ei- 
ther Arimoto-Yoshii or Salpeter), and we keep record of 
the age and the metallicity of each of these SSPs. It is 
quite straightforward to compute the global luminosities and 
colours of our simulated galaxies, as the sum of the contri- 
bution of their constituent star particles. SSP luminosities 
are computed by mass-weighted integration of the Padova 
isochrones (Girardi et al., 2002). 

We must however pay special attention to the rela- 
tion between star particle mass and SSP mass. In our "sta- 
tistical" implementation of chemical evolution (Lia et al. 
2002a,b), a fraction of the star particles with age trans- 
forms back into gas particles (the gas-again particles of Lia 
et al.), to simulate the re-ejection of gas into the interstellar 
medium by dying stars. Following the notation by Lia et al., 
the re-ejected fraction increases with the age t of the SSP: 

where $(M) is the IMF, M u is its upper mass limit, and 
M(t) is the mass of the star with lifetime t. Typical values of 
the global returned fraction after a Hubble time are 30% for 
the Salpeter IMF, and around 50% for the Arimoto-Yoshii 
IMF. As a consequence, out of an episode of star formation 
involving N star particles, after a time t on average only 
N x (1 — E{t)) remain, while N x E(t) have returned to be 
SPH particles. We need to take this effect into account when 
computing the luminosities, since SSP luminosities refer to 
the initial mass of the SSP, namely the mass involved in the 
original star formation episode, not to its current mass which 
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Table 1. Basic physical properties and prescriptions adopted for the runs of the "Virgo" and "Coma" clusters discussed in this paper. 
Temperatures (4 th column) refer to the "standard" run AY-SW. The last two columns give quantities as resulted from the simulations, 
at -Rsoo- 



cluster 


M vir 

[10 14 Mq] 


[Mpc] 


<C kT ew > 
[keV] 


IMF 


SN feedback 
efficiency 


pre-heat. @ z=3 
[keV/part.] 


therm, cond. 
(Spitzer) 


J cold 

(Rsoo) 


Z Fc 
Z Fe,Q 


"Virgo" 


2.77 


1.62 


3.0 


AY 


SW 






0.19 


0.21 










AY 


SW 




1/3 


0.17 


0.38 










AY 


SWx 2 






0.11 


0.22 










AY 


SWx 4 






0.06 


0.16 










AY 


SW 


0.75 




0.16 


0.22 










AY 


SW 


1.5 




0.15 


0.30 










Sal. 


SW 






0.23 


0.14 










Sal. 


Weak 






0.27 


0.12 


"Coma" 


12.40 


2.90 


6.0 


AY 


SW 






0.20 


0.20 










AY 


SW 




1/3 


0.17 


0.20 



is a fraction 1 — E(t) of the initial one. Each star particle 
of age t is effectively representative of 1 _ 1 E ^ star particles 
at its birth — or equivalently, each star particle of mass m* 
corresponds to an initial SSP mass j-r^m ■ Therefore, for the 
computation of the luminosity each star particle is assigned 
a corresponding "initial SSP mass" j-r^tj' ra ther than its 
actual present mass m*. 



3 STAR FORMATION HISTORIES OF THE 
CLUSTER GALAXIES 

The star formation process in the cluster environment is 
known to peak at higher redshifts than in the field. The 
morphology-density relation (e.g. Dressier 1980; Goto et al. 
2003) implies that the cluster population is dominated by 
early-type galaxies, ellipticals and SO's, whose stellar pop- 
ulations have formed at redshift z > 2 as indicated by the 
tightness and the redshift evolution of the colour-magnitude 
relation and of the fundamental plane (e.g. Bower, Lucey & 
Ellis 1992; Kodama & Arimoto 1997; j0rgensen et al. 1999; 
van Dokkum & Stanford 2003). Conversely, the field galaxy 
population is dominated by late Hubble types, with star for- 
mation histories still presently on-going. Furthermore, com- 
pared to their cluster counterparts, field ellipticals (espe- 
cially low mass ones) display more extended star formation 
histories and/or minor star formation episodes at low red- 
shifts (e.g. Bressan et al. 1996; Trager et al. 2000; Bernardi 
et al. 2003; Treu et al. 2005). The rapid drop in the past 
star formation rate in clusters is also apparent in the evo- 
lution of the fraction of star forming galaxies, as traced by 
blue colours and spectroscopy (e.g. Butcher & Oemler 1978, 
1984; Ellingson et al. 2001; Margoniner et al. 2001; Poggianti 
et al. 1999, 2004; Gomez et al. 2003). 

These trends are qualitatively in line with present the- 
ories of cosmological structure formation, predicting that 
high density regions evolve faster than low-density regions 1 

1 It is worth mentioning that even when modelling isolated 
galaxy evolution (in the monolithic collapse fashion), a more pro- 
longued star formation history is predicted in low-mass, low den- 
sity galaxies, with respect to the more massive or denser ones 
(Carraro et al. 2001; Chiosi k, Carraro 2002). 



(Bower 1991; Diaferio et al. 2001; Benson et al. 2001); 
although there are severe quantitative discrepancies with 
the results of semi-analytical hierarchical models, especially 
concerning the formation and number density evolution of 
field early-type galaxies (van Dokkum et al. 2001; Benson 
et al. 2002). 

All in all, clusters display a drastic decay in the star for- 
mation rate at z < 1, much faster than the corresponding de- 
cline indicated by the Madau plot for field galaxies (Kodama 
& Bower 2001; Fig. 1). Such drop in the cluster star forma- 
tion rate is usually ascribed to a combination of three ef- 
fects. (1) Interaction with the cluster environment and with 
the ICM quenches the star formation of the infalling galax- 
ies via mergings and interactions, harassment, ram pressure 
stripping, starvation/strangulation, etc. (Dressier 2004 and 
references therein). (2) The rate of accretion of star forming 
field galaxies onto clusters decreases at decreasing redshift 
(Bower 1991; Kauffmann 1995). (3) The intrinsic star for- 
mation rate of the accreted galaxies decreases at z < 1 (cf. 
the Madau plot, Madau et al. 1998). 

In this section, we analyze the global star formation his- 
tory of the cluster galaxies in our simulations and compare 
it to the observational one derived by Kodama & Bower 
(2001). Just like for the empirical one, our global star for- 
mation history is obtained as the sum of the star forma- 
tion histories of the individual cluster galaxies (excluding 
the young cD stars in the central regions, formed out of the 
cooling flow; see also Paper III). 

For each cluster, we single out all the star particles that 
at the present time (z=0) belong to the identified galaxies 
within 1 Mpc, corresponding to the same Rzo radius consid- 
ered by Kodama & Bower. For each star particle we know 
the age t and we derive the corresponding "initial SSP mass" 
involved in star formation at its birth time, as x™E{t) ( see 
§ 12.21 . Summing over all the selected star particles, we re- 
construct the past SFR in the comoving Lagrangian volume. 

In Fig. 1 the global normalized SFR is shown for the 
galaxies in "Virgo" (3 models with varying feedback strength 
and Arimoto- Yoshii IMF, plus a model with Salpeter IMF) 
and "Coma" (standard reference prescription). Since z ~ 1.2 
the SFR declines considerably more rapidly than that of 
the field galaxies, in line with observational reconstructions 
indicating a peak at z ~ 3 — 4 and a subsequent significative 
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Figure 1. Normalized global star formation rates in the "Virgo" 
and "Coma" galaxies, excluding the youngest stars (less than 
1 Gyr old) in the inner 10 kpc of the cDs; symbols are explained 
in section 2. Also shown are observed rates for field galaxies (dot- 
ted, Madau et al., 1998) and those of galaxies in the inner parts 
of rich clusters (shaded region, from Kodama & Bower, 2001). 

drop from z ~ 2 to quite more pronounced than in the 
field. In Fig.0 the shape of the SFR history in our simulated 
clusters is in agreement with the estimates of Kodama & 
Bower, however the normalization is somewhat higher (close 
to the field level at redshift 2 = 2). 

From comparing "Virgo" and "Coma" in the figure, it 
is evident that the total star formation history normalized 
by cluster mass is nearly independent of the cluster virial 
mass, as also recently pointed out by Goto (2005) over a 
sample of 115 nearby clusters selected from the SDSS: this 
suggests that physical mechanisms depending on virial mass 
(such as ram-pressure stripping) are not exclusively driving 
galaxy evolution within clusters. 



4 SCALING PROPERTIES OF THE CLUSTER 
LIGHT 

The near-infrared luminosity of galaxies is only negligibly af- 
fected by recent star formation activity, thus giving a robust 
measure of the actual stellar content of a cluster (Kauffmann 
& Chariot 1998). In Fig-.^we show the total (2.2fi) K-band 
luminosity of galaxies within rsoo versus the total mass in- 
side of rsoo (which in the ACDM model is approximately half 
of the virial radius) for individual clusters in the simulated 
sample — including objects of lower mass than "Virgo" and 
"Coma". The various models match quite well the best-fit 
slope derived observationally by Lin, Mohr & Stanford 2004: 
L|q oc M^qq 2 , but the runs with additional pre-heating and 
those with stronger feedback result in a normalization too 
low with respect to the observations. As to the Salpeter sim- 
ulations, the excellent agreement seen in this plot is partly 
affected by the very blue colours of the galaxies (Fig. llOU : in 
bluer bands, e.g. the B band, these simulations results too 
bright for their cluster mass (Paper I). 

The figure's insert shows the result of correcting for the 



M 500 [10» MJ 

Figure 2. .ft"- band luminosity- mass relation inside of -R500 at 2=0 
(excluding cD stars in the innermost 10 kpc, younger than 1 Gyr); 
notice that "Coma" cluster results at M500 ~ 7 X 10 14 Mq are 
available only for the AY-SW prescription, with and without ther- 
mal conduction (see Table 1). Shaded strip: best-fit relation for a 
sample of 93 clusters and groups from 2MASS (Lin et al., 2004). 
In the insert: the same excluding all stars born since z ~ 1 from 
the central 40 kpc; also the observational best-fit relation has 
been modified to exclude the BCG. 



excess of young cD stars by neglecting the luminosity con- 
tribution of stars in the innermost 40 kpc formed at z < 1 
(see Paper III for details) ; the observational data have been 
corrected as well by excluding the cluster brightest (usu- 
ally central) galaxy (BCG). The slope of the relation re- 
mains essentially the same, though observationally it steep- 
ens slightly to L|o oc M^ 8 3 (Lin et al. 2004). This indicates 
that the relative contribution of the cD/BCG to the total 
cluster light decreases with cluster mass. 

The regular trend seen in Fig. |5] over a considerable 
range of mass suggests that star/galaxy formation in groups 
and clusters tends to be approximately a scaling process; 
moreover the slope of the observational relation indicates 
that smaller groups produce stars at higher efficiency than 
larger clusters, provided that K-band light well traces stellar 
mass. 



5 THE LUMINOSITY FUNCTION OF 
CLUSTER GALAXIES 

In this section we discuss the luminosity function of the pop- 
ulation of galaxies in our simulated "Virgo" and "Coma" 
clusters (excluding the central cDs, which are discussed 
in Paper III). As mentioned above, we identify for the 
standard runs 42 galaxies in "Virgo" and 212 galaxies in 
"Coma" , respectively, with a resolution-limited stellar mass 
> 2.5 x 10 9 Mq, of the order that of the Large Magellanic 
Cloud. This limit corresponds to Mb ~ —16 at z=0 (Fig. 3). 

In Fig. |3]we compare the B-band and K-band luminos- 
ity function (LF) of our simulated cluster galaxies, to the 
observational LFs by Trentham (1998; a weighted mean of 9 
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Figure 3. Top panels: B-band luminosity functions at 2=0 for 
the standard (AY-SW) "Virgo" (left) and "Coma" (right) simula- 
tions; the LFs are normalized to have the same number of galaxies 
as the observed average LF (Trentham 1998, triangles) within the 
populated luminosity range (—19 $C Mg $C —17). Bottom panels: 
K-band LFs for the same clusters, compared to the observed av- 
erage LF by Lin et al. (2004); the dotted line marks the limit 
Mk ^ —21 where the observational estimate is considered reli- 
able. 



low-redshift clusters), and by Lin et al. (2004; a composite 
LF of 93 clusters and groups from 2MASS data). 

For the B-band LF (top panels), the number of resolved 
individual galaxies quickly drops for objects fainter than 
Mb ~ — 17, both in "Virgo" and in "Coma", due to the 
above mentioned resolution limits and identification proce- 
dure. Notice that, although we are missing the dwarf galax- 
ies that largely dominate in number, we are able to describe 
the bulk of the stellar mass and of the luminosity in clus- 
ters, which is dominated by galaxies around L* while dwarfs 
fainter than Mb ——17 give a negligible contribution (Cross 
et al. 2001; Blanton et al. 2004) — though in the simula- 
tions, too much of the luminosity and stellar mass expected 
in objects around L t , is instead accumulated in the cen- 
tral cD (see below). Our B-band LF in Fig. [3] is normal- 
ized so that, in the luminosity bins that are significantly 
populated ( — 19 ^ Mb ^ — 17), the overall number of ob- 
jects is the same as for the observed relative LF (Trentham 
1998); in this magnitude range, the shape of the predicted 
and observed LF is directly comparable. The simulated LF is 
steeper than the observed one, namely we underestimate the 
relative number of bright galaxies. One reason for this may 
be our somewhat biased selection of the group and clusters 
sample, which consists of fairly relaxed systems, in which no 
significant merging takes place since z < 1 (see Paper III). 
This means that most of the massive galaxies have already 
merged with the central cD by dynamical friction - which 
is not always the case for real clusters (one likely, nearby 
example of an unrelaxed cluster is Virgo, e.g., Binggeli et 
al. 1987). 

In the K-band (bottom panels of Fig.[3Jl, our simulated 
LF reaches fainter magnitudes than the limit Mk < —21 
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Figure 4. Same as in top panels of Fig. [3] for various "Virgo" 
simulations; symbols are explained in section 2. 



(dotted line), where the observational LF is considered re- 
liable (Lin et al. 2004). Therefore, we normalize our LF to 
match the observed number of galaxies within the populated 
magnitude range and above the resolution limit (—23.2 ^ 
Mb ^ —21). In the K-band, the lack of massive galaxies in 
relative number is even more evident. 

In Fig. [I] we compare the LFs for the "Virgo" clus- 
ter simulated with different IMFs and feedback prescrip- 
tions. For the three AY simulations, increasing the feedback 
strength results in decreasing the masses of all galaxies, and 
hence, in particular, the number of bright galaxies. The re- 
sults with thermal conduction (not shown in the figure) are 
very close to the analogous case without thermal conduction. 
The LF in the Salpeter case is broader in luminosity, due the 
lower stellar feedback with respect to the AY IMF, allow- 
ing a larger accumulation of stellar mass in galaxies. The 
trend is even stronger in the Sal-WFB case. Notice that, 
though the Salpeter IMF simulations are more successful 
in forming massive galaxies and would thus compare more 
favourably to the observed LF, this comes at the expense of 
a too large cold fraction with respect to observations (Ta- 
bleland Paper I). Other reasons to disfavour the Salpeter 
IMF simulations are the low predicted metallicities in the 
ICM (Tableland Paper I), and the too blue colours of the 
galaxies (§EJ. 

In Fig. |S]we compare the LF of the "mini-Coma" clus- 
ter, in absolute number of galaxies, to the observed Coma 
LFs: besides discussing the distribution of the stellar mass 
in the simulations (cD vs. bright galaxies vs. dwarf galax- 
ies) we want to test if the actual number of the galaxies 
formed in the simulation is sensible. The Coma LF is from 
Andreon & Cuillandre (2002) and Andreon & Pello (2000); 
for the comparison, we compute the luminosity functions of 
galaxies in our simulated cluster within projected, off-center 
areas comparable to the areas covered by the observations 
- about (0.8 Mpc) 2 in B, (1 Mpc) 2 in V, R and (0.5 Mpc) 2 
in H. There is quite good agreement with the observed num- 
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Figure 5. Absolute luminosity functions at z=0 of the "Coma" 
cluster within the relevant projected area (see text) in the B, V, 
R, H bands, compared with observed Coma LFs from Andreon 
& Cuillandre 2002 (B, V, R), and Andreon & Pello 2000 (H). 

ber of galaxies over the luminosity range that is significantly 
populated by the galaxies identified in the simulation. How- 
ever, as already noticed from the relative LF in Fig. |3] we 
are "missing" the bright end tail of the LF, a problem which 
does not seem to be cured by increased resolution (see H5.ip . 
These results on the absolute LF confirm that the prob- 
lem is not that we form too few galaxies in the simulations 
(whence one could statistically "miss" the fewer galaxies at 
the bright end of the LF), but that in fact there is a lack of 
bright galaxies in favour of the overgrown cD. 

5.1 Resolution effects 

In order to test for numerical resolution effects, one simula- 
tion of the "Virgo" system with the standard AY-SW model 
was run at 8 times higher mass and 2 times higher force res- 
olution. At the time of writing, this simulation has reached 
z ~0.8, so that in Fig. we can compare the Virgo LF 
for the two different resolutions at z ~0.8. As expected, at 
higher resolution the LF extends to fainter magnitudes. This 
is due to the higher force and mass resolution, which allows 
one to resolve as several small (proto-)galaxies an object 
which was previously initially identified as one large proto- 
galaxy. However, above the resolution limit for the galaxy 
identification, the LF appears to be strikingly robust to res- 
olution effects, and consequently so is all the discussion in 
the previous section. The global stellar content of the clus- 
ter is also essentially resolution independent: the total mass 
of stars at z ~0.8 inside of the virial radius is 3.0 ■ 1O 12 M0 
for the normal resolution case and 3.4- 1O 12 M for the high- 
resolution one. The corresponding masses of stars in galaxies 
other than the cD are 9.1 and 9.2-10 11 Mq; the total masses 



of stars outside of the cD (adopting r C D=50 kpc) are 1.7 
and 2.OTO 12 M0, and the masses of stars in the cD are 1.3 
and 1.4-10 12 M©, respectively. The total numbers of galaxies 
identified, assuming N m in ~7 in both cases, are 50 and 341, 
respectively. 

In conclusion, at higher resolution more substructure 
can be identified extending the LF at the dwarf galaxy end; 
yet the luminosity of the second brightest galaxy (after the 
cD) in the two simulations of different resolution is quite sim- 
ilar. Moreover, the number of galaxies brighter than « L* 
remains fairly unaffected when going to higher resolution: 
enhancing resolution does not increase the number of bright 
cluster galaxies. Analogous results are obtained from an- 
other resolution test we performed, on a smaller system (a 
group of galaxies of virial mass 0.48 -1O 14 M ) evolved down 
to z = 0: the higher resolution does not affect the luminosity 
of the second ranked galaxy nor helps increasing the number 
of bright (~ L, ) galaxies. These tests indicate that the lack 
of bright galaxies in the normal resolution simulations of our 
"Virgo" and "Coma" clusters is not an effect of numerical 
"over-merging" . 

"Classic" semi-analytical models of cluster formation 
also faced, to some extent, the problem that the cD became 
too bright and the number of lower ranked, bright galaxies 
was depleted. Springel et al. (2001) suggested that numerical 
over-merging onto the central cD in dissipationless dark 
matter simulations could deplete the number of galax- 
ies at the bright end of the LF, and showed that for semi- 
analytical models this problem can be significantly improved 
on by suitable identification of dark matter substructure, or 
of "haloes within haloes". This is not a solution in our case, 
however, because we identify as galaxies all bound systems 
containing as little as jV min =7 star particles, and increasing 
resolution does not increase the number of bright galaxies. 
Moreover, in hydrodynamical simulations there is as of yet 
no physical way to quench the central, semi-steady cooling 
flow; such late gas accretion accounts for part of the excess 
of the cD masses. 
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5.2 Redshift evolution of the Luminosity Function 

The LF in clusters can evolve due to two effects: passive lu- 
minosity evolution from the aging of the stellar populations, 
and mass evolution due to dynamical effects such as mergers, 
tidal stripping and dynamical friction. Cluster galaxies are 
known to consist of old stellar populations, with the bulk of 
their star formation occurring at z^>l. However, they may 
not have been completely assembled by 2 ~ 1, and they 
may still undergo one or two significant mergers since then, 
or grow by gas accretion (van Dokkum et al., 1999). 

The observed evolution of the luminosity function is a 
powerful probe for the assembly epoch of galaxies. From 
a theoretical point of view, hierarchical models of galaxy 
formation based on semi-analytical prescriptions tend to 
predict a deficit of massive galaxies at z ~ 1, as a re- 
sult of the ongoing mass assembly activity at lower redshift 
(Kauffmann & Chariot, 1998). More recent models, based 
on a flat ACDM cosmology and with updated star forma- 
tion, feedback and wind prescriptions perform much better 
and agree with the observed number of massive galaxies up 
to z — 1.2 — 1.5; however they still face problems with a 
deficit at higher redshifts as well as with the colours of mas- 
sive galaxies, underpredicting the number of red early type 
galaxies and of Extremely Red Objects. In the deep field, the 
observed evolution seems to lie in between the predictions 
of current hierarchical models and those of the monolithic, 
pure luminosity evolution scenario (Somerville et al. 2004; 
Stanford et al. 2004). 

For clusters instead, recent results from the observa- 
tions of distant objects suggest that the evolution of the 
LF is best modelled by pure luminosity evolution: Barger et 
al. (1996, 1998), De Propris et al. (1999), Kodama & Bower 
(2003) found that, after correcting for this effect and for cos- 
mological dimming, the galaxy mass function has actually 
evolved little over time; the K-band LF at 2 ~ 1 is con- 
sistent with pure luminosity evolution with constant stellar 
mass and early redshifts of star formation (zf > 2). 

In Fig. |7| we show the evolution of the B and K band 
LFs of our "Virgo" and "Coma" clusters. We computed the 
(rest-frame) LF of the galaxies identified in the z = 1 frame; 
this is shifted to brighter magnitudes with respect to the 
distribution at z = 0, and both luminosity evolution and 
dynamical mass evolution are seen to play a role in our sim- 
ulated clusters. In particular, since the bulk of the stars in 
our galaxies are formed at z > 2 (see Paper III and Fig. |3] 
bottom panel), luminosity dimming is an important effect. 

We also computed the LF at z — 1 as expected from 
pure passive evolution of the z — LF, by considering the 
star particles in each of the galaxies identified at z = 0, 
and computing the luminosity they had back at a ~7.5 Gyr 
younger age. Such "pure passive evolution" LF is displayed 
as a dotted line in Fig. [7] and it matches quite closely the 
actual LF in the z = 1 frame (apart from the absence of a few 
bright galaxies, see below). This indicates that luminosity 
dimming of the stellar populations drives most of the LF 
evolution and the fading to fainter magnitudes from z — 1 
to z — for the simulated galaxy population; this is in line 
with the results of Kodama & Bower (2003). 

In the brightest luminosity bins of our LF, some ad- 
ditional, dynamical effect seems to be required since more 
bright galaxies are identified at z = 1 than expected from 
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Figure 7. Top panels: Luminosity functions in absolute number 
of galaxies per B-band magnitude bin for the standard (AY-SW) 
"Virgo" and "Coma" simulations, at 2 = (solid lines) and 2 = 1 
(dashed lines); the dotted lines represents the expected z = 1 LF 
if pure passive evolution of the stellar populations is applied to 
the z = LF. Middle panels: Same, for the K-band luminosity 
functions. Bottom panels: Mass function of cluster galaxies at 
2 = and 2 = 1. 



pure passive evolution; especially in the case of the "Coma" 
simulation, where 263 galaxies ar found at z — 1 versus 212 
at 2 = 0. To assess mass evolution effects, we plot in the 
lower panels of Fig.|7|the mass function of (the stellar com- 
ponent of) cluster galaxies at 2 = 1 and at 2 = 0. The two 
mass functions are very similar at the low mass end, but 
a few rather massive objects (above M* ~ 4 x 10 10 Mq) 
present at 2 = 1 "disappear" at 2 = 0. Inspection of the 
simulations show that these objects have merged with the 
cD at 2=0. 

As shown in § 15.11 numerical "over-merging" onto the 
central cD cannot be the main reason for missing galaxies at 
the bright end of the simulated LF. The lack of bright, mas- 
sive galaxies other than the cD at 2=0 is probably largely 
affected by our group and cluster selection procedure, which 
picks out extremely relaxed objects at z=0, as discussed in 
SjK] All in all, "Coma" presents a less relaxed structure than 
"Virgo", as also inferred by the shapes of the LF at z=0 
and z=l. Indeed even for "Virgo" we do see a deficiency 
of bright galaxies at z=0, but not at z=l, where a gap is 
just due to poorer statistics of galaxy numbers. Work is in 
progress however to analyze a large sample of groups more 
randomly selected (hence including both relaxed and un- 
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relaxed object) to assess the bias induced by the selection 
procedure (Sommer-Larsen, D'Onghia & Romeo 2005, in 
preparation). 

5.3 Passive evolution and the Fundamental Plane 

In the previous section we computed the pure passive evo- 
lution dimming, from z = 1 to z = 0, corresponding to the 
cluster galaxies identified at z — 0. Constraints on the pas- 
sive evolution of ellipticals in clusters are derived from the 
observed evolution of the Fundamental Plane (FP) , indicat- 
ing a dimming of about 1.2 magnitudes in the B band in the 
redshift range z = 1 — > 0; assuming a Salpeter slope for the 
IMF, this implies a redshift of formation zf or ^ 2.5 for the 
stellar populations (Van Dokkum & Stanford 2003; Wuyts 
et al. 2004; Renzini 2005; Holden et al. 2005). 

A direct comparison to FP constraints is hampered by 
the fact that in our LF we miss the bright, massive elliptical 
galaxies defining the observed FP; nevertheless, it is inter- 
esting to comment on the predicted passive evolution of our 
galaxies. In the table below we show the B magnitude evo- 
lution AMb = Mb(z = 0) - M B (z = 1) predicted by the 
SSPs in use, as a function of the assumed redshift of forma- 
tion Zf or of the stellar population. These values are largely 
independent of the SSP metallicity, at least within a factor 
of 3-4 of the solar value (which is certainly representative 
of the bright galaxies in our simulations) . 

zf or AMb (Salp) AM B (AY) 



5 
3 

2.5 
1.5 



0.97 
1.15 
1.26 
2.01 



1.07 
1.27 
1.40 
2.20 



Since younger stellar populations dim faster, the lower the 
redshift of formation, the faster the magnitude evolution be- 
tween z — 1 and z — 0. For the Salpeter IMF, our photomet- 
ric code predicts that the observed dimming of ~1.2 mag, 
corresponds to Zf or ^2.5 in agreement with the above men- 
tioned studies. The rate of dimming depends however also 
on the slope of the IMF, being faster for shallower IMFs; 
for an AY slope, the luminosity evolution is in fact slightly 
faster requiring Zf OT ^ 3 (see also Renzini 2005). Taking 
into account progenitor bias, i.e. the fact the progenitors 
of the youngest present-day early-type galaxies drop out of 
the sample at high redshift, the actual magnitude evolution 
of ellipticals might be underestimated by about 20% (van 
Dokkum & Franx 2001). In this case, AMb = 1.4 mag and 
Zf or = 2.5 (AY IMF) are still compatible with observations. 

The mean redshift of formation for the stars in our simu- 
lated galaxies is 2.5 (Paper III), which is on average compat- 
ible with FP constraints but the scatter around the mean is 
large. Younger star particles, where present, dominate in the 
luminosity contribution at z = 1 and induce a much faster 
overall evolution (cf. a stellar population formed at z=1.5 
dims by more than 2 mag). This scatter in age makes some 
of our galaxies evolve faster than indicated by FP studies. 

This is illustrated in Fig. |H| where we show the mag- 
nitude evolution AMb = Mb[z — 0) — Mb (z = 1) for our 
simulated "Coma" galaxies as a function of their present- 
day magnitude and as a function of average stellar age. The 
age-dimming (i.e., Zf OT — AMb) relation is well defined in 
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Figure 8. Passive evolution in B magnitude between z = 1 and 
z = 0, predicted for the Coma cluster galaxies. AMb = Mb(z = 
0) — Mb(z = 1) is plotted as a function of the present— day mag- 
nitudes at z = (top panel) and as a function of age (bottom 
panel). For the six brightest objects (largw symbols) we also show 
the internal spread in stellar ages responsible for the very large 
magnitude evolution (see text). 



the bottom panel and the FP constraint of 1.2-1.4 mag is 
met by objects older than 10.5-11 Gyr. The scatter above the 
lower envelope of the age-dimming relation is due to inter- 
nal age scatter within the galaxies. We analyze in particular 
the six brightest objects in the simulation, with Mb < — 19 
(large symbols in the figure), which are more relevant for 
comparison to the bright spheroidals on the FP. These ob- 
jects are 9-10 Gyr old (average age), i.e. z/ or =1.5-2 which 
corresponds to AMb = 1.5 — 2. However, they also present 
a large internal scatter in stellar age, as shown by the er- 
rorbars in the bottom panel (stretching from the minimum 
to the maximum stellar age in each object). In particular, 
they contain stars as young as 8 Gyr, and some of them have 
minor tails of SF stretching to z < 1, i.e. to ages younger 
than 7.5 Gyr. The presence of these younger-than-average 
stars has negligible effects on the present-day appearance of 
these galaxies, since at z = all the stars are already quite 
"old" (ages > 5 Gyr) — see also Fig.|UJI showing that the lu- 
minosity age equals the actual average age for these objects, 
while it would be significantly younger if the younger stellar 
component were prominent. However, back at z = 1 these 
younger-than-average stars become very bright, young pop- 
ulations (ages ^1 Gyr) with a major contribution to lumi- 
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nosity, so that they ultimately drive the overall brightening 
of the galaxy when we trace it back to to z = 1. 

The median of the luminosity evolution for all the Coma 
galaxy population is 1.4 mag (in agreement with the average 
zf or = 2.5), i.e. close to the FP constraint; however the 
large internal age scatter within the galaxies, especially the 
brighter ones, induces an overall predicted passive evolution 
for the LF more extreme than that (Fig. |TJ| . 

Higher resolution simulations can possibly help with 
this problem, by both anticipating the overall SF process (re- 
solving smaller, hence denser, substructures) and reducing 
the internal poissonian noise in the stellar age distribution 
of the individual galaxies. However, the high resolution test 
presented in § l5.1l does not show major differences in the LF 
down to z = 0.8 (Fig. EJ, suggesting that we have reached a 
resolution good enough to grasp the main galaxy properties. 
We are probably then facing here a standard problem of hi- 
erarchical models of galaxy formation, i.e. that they predict 
both slightly younger average ages (Fig. |3Jl) and a larger 
(internal) age scatter for more massive objects, both trends 
in contrast with observational evidence (see Section [3] and 
references therein). 

We notice however that the faint end of the present-day 
colour-magnitude relation of E+SO cluster galaxies seems 
to be largely populated by objects which reddened onto the 
Red Sequence only recently, while at high z they were much 
bluer; an effect that is usually associated to the spiral — > 
SO transformation and to the Butcher-Oemler effect from 
intermediate redshift to z = (e.g. De Lucia et al. 2004; 
see also Section[3]and references therein). These objects are 
certainly not the same probed by the FP at high redshift, 
which relates to the most massive ellipticals. Also, Holden 
et al. (2005) find a large scatter in mass-to-light ratio and 
hence magnitude evolution for the least massive ellipticals 
selected at high z, with masses around 10 11 Mq (comparable 
to the most massive objects in our final galaxy population, 
Fig. . Henceforth, it is possible that our simulations miss 
in fact exactly the massive spheroidals at the bright end of 
the LF, which are most meaningful for FP constraints. The 
issue needs then to be revisited once the problem of properly 
populating the bright end of the LF is solved for simulated 
clusters. 



6 THE RED SEQUENCE 

The light and stellar mass in clusters of galaxies is domi- 
nated by bright, massive ellipticals (Abell 1962, 1965). These 
are known to form a tight colour-magnitude relation, or 
Red Sequence (Bower, Lucey & Ellis 1992ab; Gladders el al. 
1998; Andreon 2003; Hogg et al. 2004; Mcintosh el al. 2005). 
In Fig. EJib we compare the colour-magnitude relation for 
our "Coma" cluster galaxies, excluding the cD, to the ob- 
served Red Sequence of Coma from Bower el al. (1992) and 
Terlevich, Caldwell & Bower (2001), assuming a distance 
modulus to Coma of 35.1, i.e. Ho ~70 km sec -1 Mpc -1 
(Baum el al. 1997). The slope of the observed Red Se- 
quence is sensitive to aperture effects: fixed-aperture obser- 
vations are often used, measuring a smaller fraction of the 
total light in larger objects than in smaller ones; combined 
with colours gradients, this introduces a bias in the sense of 
measuring systematically redder colours in larger galaxies, 



i.e. a steepening of the actual colour-magnitude relation. 
When the analysis is instead performed on an area that 
scales with the size of the galaxy, e.g. within the effective 
radius, or within a given isophotal/photometric radius, the 
colour-magnitude relation is in fact flatter (Terlevich el al. 
2001; Scodeggio 2001; and references therein). For Coma and 
Virgo, Bower et al. measured colours within a fixed aperture 
of 11" and 60" respectively, corresponding to a fixed radius 
of 5/i -1 ~ 7.15 kpc. Henceforth for a proper comparison, we 
computed and plotted in Fig.|5|the colours for our simulated 
galaxies within a projected radius of 7.15 kpc. The magni- 
tude in abscissa are instead total magnitudes, as adopted by 
Bower et al. 

The solid lines in panels (a,b) indicate the observed 
slope and location of the Red Sequence; the dotted line is a 
least square fit to our data. There is overall good agree- 
ment, although our Red Sequence is not as extended as 
the observed one (which reaches magnitudes as bright as 
Mv = —23), due to the lack of simulated galaxies at the 
bright end of the luminosity function discussed in the previ- 
ous section. Our slope for the (U-V) vs. V relation (Fig.|UJi), 
is slightly flatter than the observed one, but the colours of 
the three brightest objects are in good agreement with ob- 
servations, even within the observed small scatter. These 
objects are the most significant ones if we remember that 
the observed slope is defined only for My < —19, extended 
down to Mv < —18 by Terlevich et al. The slope of the sim- 
ulated (V-K) vs. V relation (Fig.^) is in perfect agreement 
with the observed one. 

Although we formally obtain the right slopes, the simu- 
lated Red Sequence also displays a much larger scatter than 
the observed Coma Red Sequence, which is very narrow 
down to Mv ~ —18 (dashed and dotted lines in Fig. |5^b; 
see also Fig. 2a in Gladders et al., 1998). Such large scatter 
is partly due to Poissonian noise from the small number of 
star particles in the low luminosity galaxies. Besides, while 
the observed colour-magnitude relation refers to early type 
galaxies (E/S0), we did not attempt any morphological clas- 
sification and a priori selection of our simulated galaxies due 
to limited resolution. 

The colour-magnitude relation is classically interpreted 
as a mass-metallicity relation (Kodama & Arimoto 1997). 
In our simulations, the bulk of the stars in galaxies are 
formed at z >2 (see also Paper III), hence age trends 
among galaxies are only mild (Fig. |3Jl) and the colour- 
magnitude relation is mostly a metallicity effect (Fig. ^J;). 
The metallicity-luminosity relation for our cluster galax- 
ies is displayed Fig. |3J;, both in terms of mass-weighted 
metallicity and of (V-band) luminosity-weighted metallic- 
ity. The latter is systematically slightly lower than the actual 
mass-averaged stellar metallicity, as expected since more 
metal-poor populations tend to be brighter, skewing the 
luminosity-weighted metallicity to lower values (Greggio 
1997); the trend with galaxy luminosity is however main- 
tained. Although metal-rich objects seem to exist at all lu- 
minosities, the fraction of metal-poor galaxies (as well as the 
scatter in metallicity) increases with decreasing galaxy lu- 
minosity; the average stellar metallicity decreases for fainter 
galaxies (dotted line), which ultimately drives the simulated 
colour-magnitude relation. 

Fig. |nji shows the average stellar ages for the galax- 
ies. All the galaxies consist of old stellar populations, with 
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Figure 9. (a— b) Colour-magnitude relation for galaxies in the 
simulated "Coma" cluster (data points and dotted linear fit line), 
compared to the observed relation with its scatter (solid and 
long-dashed lines, from Bower et al. 1992; short-dashed lines, 
from Terlevich et al. 2001). The extension of the dashed lines 
indicates the magnitude range probed obscrvationally, while the 
average fit relation (solid line) has been extrapolated to lower 
magnitudes, (c) Metallicity— luminosity relation for the "Coma" 
galaxies (full symbols for mass-averaged stellar metallicity, open 
symbols for luminosity-weighted metallicity); the dotted line is a 
linear fit. (d): Age— luminosity relation for the "Coma" galaxies 
(mass-average and luminosity— weighted stellar ages); the dotted 
line is a linear fit. 



average ages between 7 and 12 Gyrs. The brighter galax- 
ies (My < —20) are essentially coeval, 9-10 Gyr old. The 
luminosity-weighted ages are generally very close to the 
mass averaged ages, which implies a small scatter of stel- 
lar ages within the individual galaxies: extended tails or re- 
cent episodes of star formation, would strongly skew the 
luminosity-weighted estimate to younger ages. The scatter 
in age among the galaxies is around 30%, with a mild trend 
of decreasing age with increasing luminosity (dotted line), 
but the effect is much smaller that the systematic variation 
in metallicity (a factor of 10, Fig.^). This confirms that the 
color-magnitude relation of our simulated galaxies is driven 
by a mass-metallicity relation, in agreement with common 
wisdom. The mild age trend, which tends to act in the op- 
posite direction (making fainter galaxies redder) is probably 
responsible of the fact that our colour-magnitude relation is 
not as steep as the observed one. 

At all magnitudes, even the faintest ones, we find some 



simulated galaxies with large (super-solar) metallicities. For 
dwarf ellipticals, a very large spread in metallicities and 
colours with a tail of red, metal rich objects, is in fact ob- 
served but at fainter magnitudes than those probed here 
[M B > -15, Conselice et al. 2003). Metal rich dwarf el- 
lipticals in clusters may possibly originate from larger (and 
hence more metal rich) galaxies that have suffered consid- 
erable tidal stripping in the cluster potential. To test this 
hypothesis we traced the evolution of the nine galaxies with 
Z, > Zq and My > —18, shown in the top panel of Fig. 8, 
back in time. Two of the galaxies were found to contain a 
significantly larger mass of stars at z ~ 1.5 — 2.5 (when 
these galaxy reach their maximum stellar masses), whereas 
the stellar masses of the other seven were found not to 
change much since their "birth". Therefore, in our simula- 
tions, stripping of originally larger galaxies is one possible, 
but not predominant, channel to form dwarf galaxies of high 
metallicity. 

In Fig. 1101 we assess the dependence of the simulated 
Red Sequence on the parameters of the simulation, by con- 
sidering the "Virgo" cluster computed with different physi- 
cal prescriptions (see Section[5]and Table0. For all the AY 
simulations the galaxies scatter about the observed Red Se- 
quence, and their average location in the colour-magnitude 
diagram appears to be robust with respect to the inclusion 
of thermal conduction or preheating, and also to enhanced 
feedback efficiency (AY-SWx2; we do not display the case 
AY-SWx4 since the total number of galaxies and their av- 
erage mass are quite small — see Fig. 0] — but they still 
scatter around the observed line). The scatter is however 
very large and the magnitude range of the galaxies formed 
quite limited (missing the bright end of the LF), so that one 
cannot make meaningful comparisons to the observed slope 
of the Red Sequence. 

Significantly different is the case of the Salpeter IMF 
simulations. The scatter is much reduced and the slope of 
the observed colour-magnitude relation is well reproduced; 
however, the simulated galaxies are offset to the blue of the 
observed Red Sequence. Let us discuss the latter problem 
first. In the standard SuperWind case (Sal-SW) in partic- 
ular, the stellar populations are too metal poor, hence too 
blue, to match the observed Red Sequence. Notice that the 
blue colours are not due to younger stellar ages, since the 
star formation history in the Virgo cluster is very similar for 
the Salpeter and the AY simulations (Fig. 0. In the weak 
feedback scenario, where effective supernova energy injec- 
tion is limited to the early epochs (Sal-WFB) the simulated 
Red Sequence falls closer to the observed one, though still 
on the blue side. With minor feedback the galaxies retain 
a much larger fraction of the metals they produce, thereby 
the stellar metallicities can grow larger; however, the price 
is a predicted enrichment of the ICM far below the observed 
levels for the Sal-WFB simulations (Tableland Paper I). 
Henceforth, top-heavier IMFs than the Salpeter IMF are 
needed, not only for the sake of the ICM enrichment but 
also to reproduce the observed colours and metallicities of 
the stellar populations. In fact the Salpeter IMF does not 
produce enough metals, and/or locks too much mass in stars, 
to account for the observed metallicities in the ICM and in 
cluster galaxies at the same time (Portinari et al. 2004, and 
Paper I). 

As to the increased dispersion in galaxy properties for 
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Figure 10. Red Sequence for our simulated "Virgo" cluster with 
different physical input. AY: Arimoto & Yoshii IMF; Sal: Salpeter 
IMF; SW: "standard" feedback prescription; SWx2: two times 
stronger feedback efficiency; WFB: weak feedback (strong feed- 
back active only at early times); COND: thermal conduction in- 
cluded; ph: energy preheating at z = 3. Symbols follow the metal- 
licity coding of Fig. [5] 



the AY vs. Salpeter simulations, this is likely induced by the 
much stronger feedback (a combination of the SW prescrip- 
tion with a top-heavy IMF). As a consequence, the Red 
Sequence in the AY simulations is much less tight, with 
no well-defined slope, although the average colours of the 
galaxies agree with observations much better than in the 
Salpeter case. If the dispersion is mostly due to numerical 
effects, because of the small number of particles in the indi- 
vidual galaxies, higher resolution simulations should show a 
reduced scatter — something to be tested in future work. 

In summary, the zero-point of the simulated Red Se- 
quence appears to be a robust prediction of the simulations, 
quite unaffected by the adopted physical prescriptions other 
than the chosen IMF, which sets the typical stellar metallic- 
ity attainable in cluster galaxies. The observed slope of the 
colour-magnitude relation is well reproduced in the Salpeter 
simulations (plagued however by an offset to too blue colours 
and low metallicities); for the AY simulations the scatter is 
much larger but the average colours are better reproduced. 



In this paper we have presented for the first time and ana- 
lyzed the properties of the galaxy population in clusters, as 
predicted from full ab initio cosmological + hydrodynamical 
simulations. 

Our results are based on cosmological simula- 
tions of galaxy clusters including self-consistently metal- 
dependent atomic radiative cooling, star formation, super- 
nova and (optionally) AGN driven galactic super-winds, 
non-instantaneous chemical evolution, effects of a meta- 
galactic, redshift dependent UV field and thermal conduc- 
tion. In relation to modelling the properties of cluster galax- 
ies this is an important step forward with respect to previ- 
ous theoretical works on the subject, e.g. based on semi- 
analytical recipes super-imposed on N-body only simula- 
tions. 

The global star formation rates of the "Virgo" and 
"Coma" cluster galaxies are found to decrease very signifi- 
cantly with time from redshift z=2 to 0, in agreement with 
what is inferred from observations of the inner parts of rich 
clusters (e.g., Kodama & Bower 2001). 

We have determined galaxy luminosity functions for the 
"Virgo" and "Coma" clusters in the B, V, R, H and K bands; 
the comparison to observed galaxy luminosity functions re- 
veals a deficiency of bright galaxies (Mb<-20). We carried 
out a test simulation of "Virgo" at eight times higher mass 
resolution and two times higher force resolution; the results 
of this test, still running at the present, indicate that the 
above mentioned deficiency of bright galaxies is not due 
to "over-merging" ; higher resolution simulations of "Coma" 
clusters are in progress as well to further check this point. 
From a suite of simulations for the "Virgo" cluster with dif- 
ferent input physics, we find that the deficiency of bright 
galaxies becomes less prominent with decreasing super-wind 
strength, in particular for models invoking a Salpeter IMF 
and only early feedback; in fact more mass can be accumu- 
lated in stars and galaxies, with low feedback. Such models, 
however, present various drawbacks: the cold fraction is too 
high and the metal production is too low, as seen in the too 
blue colours of the galaxies and/or in the low metallicity of 
the ICM, which can hardly be enriched to the observed level 
of about 1/3 of solar abundance; the latter point is discussed 
in detail in Paper I, but it also follows from more general ar- 
guments (Portinari et al. 2004). The bright galaxy deficiency 
might be explained as a selection effect, in the sense that we 
have selected cluster haloes for the TreeSPH re-simulations, 
which are "too relaxed" compared to an average cluster halo, 
so that the brightest galaxies have by now merged into the 
central cD by dynamical friction; we shall return to this in 
a forthcoming paper. 

The redshift evolution of the luminosity functions from 
redshift z—1 to is mainly driven by passive luminosity 
evolution of the stellar populations, but also by the above 
mentioned merging of bright galaxies into the cD. 

The slope of the colour-magnitude relation of the sim- 
ulated galaxies is in good agreement with the observed 
one, however the scatter is larger than observed, partly due 
to poissoinian noise within the fainter galaxies which are 
formed by small numbers of star particles. Such internal 
dispersion in stellar ages is also responsible for a luminosity 
dimming between z = 1 and z = 0, faster than indicated by 
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the observed evolution of the Fundamental Plane. The typ- 
ical average galaxy colours are best matched when adopt- 
ing a top-heavy IMF (as originally suggested by Arimoto & 
Yoshii 1987), while Salpeter IMF simulations yield too blue 
colours. Moreover we find that the average metallicity of the 
simulated galaxies increases with luminosity, and that the 
brightest galaxies are essentially coeval. Hence, the colour- 
magnitude relation results from metallicity rather than age 
effects, as concluded by Kodama & Arimoto (1997) on the 
basis of its observed evolution. 
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